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Abstract 



System of Dirac fermions with random-varying mass is studied in detail. We refor- 
mulate the system by transfer-matrix formalism. Eigenvalues and wave functions are 
obtained numerically for various configurations of random telegraphic mass m(x). 
Localized and extended states are identified. For quasi-periodic m(x), low-energy 
wave functions are also quasi-periodic and extended, though we are not imposing 
the periodic boundary condition on wave function. On increasing the randomness of 
the varying mass, states lose periodicity and most of them tend to localize. At the 
band center or the low-energy limit, there exist extended states which have more than 
one peak spatially separate with each other comparatively large distance. Numerical 
calculations of the density of states and ensemble averaged Green's functions are ex- 
plicitly given. They are in good agreement with analytical calculations by using the 
super symmetric methods and exact form of the zero-energy wave functions. 

1 Introduction 

Random disordered system is one of the most important problem in condensed matter 
physics. In two or lower dimensions, it is expected that almost all states are localized 
by the existence of random potentials, and study of these random systems obviously 
requires non-perturbative methods. In most of cases there exists no controllable pa- 
rameter for analytical calculation. However recently, exact results have been obtained 
in certain interesting disordered systems in low dimensions |T|, 0, [|. There almost all 
states are localized as expected from the general consideration, and extended states 
exist only at isolated points in physical parameter region. Mathematical tools like 
super symmetric (SUSY) methods etc. play an important role in averaging quenched 
disorders. 

The random hopping tight-binding (RHTB) model in one dimension is one of the 
most extensively studied model 0, f|, |5|. This model is related with a random spin 
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chain, a doped spin ladder, etc. ||], ||. By using SUSY methods, the density of states 
(DOS) and single-fermion Green's functions were obtained for white-noise random 
hopping||. Recently, the case of non-locally correlated disorders was studied by 
extending SUSY methods]^, §. This generalization is practically important because 
in almost all real systems random disorders have nonlocal correlations. 

In this paper we shall revisit the system of random-mass Dirac fermion, which 
describes low-energy excitations near the band center of the RHTB model. However 
our approach in this paper is somewhat different from those given so far, and we hope 
that our study and previous works are complementary with each other. In most of 
studies, average over quenched disorder is taken in order to obtain physical quantities. 
In this paper however, we shall first study quantum mechanical states in various 
fixed random disorder backgrounds. To end this, we use the transfer matrix (TM) 
formalism. The TM formalism was used to calculate the transmission and reflection 
coefficients of current incident to a one-dimensional disordered system having a certain 
momentum^]. In this paper, we shall extend the TM formalism and apply it to bound 
states and/or localized states to find their energy eigenvalues and wave functions. On 
controlling randomness, we can see how states tend to localize and what kind of states 
remain extended and contribute to long-distance correlations. To this end, results of 
our previous study on non-locally correlated disorder case is quite useful. After this 
observation, we shall take an ensemble average with respect to random variables and 
calculate the DOS and Green's functions. We show that numerical calculations are 
in good agreement with analytical results in Refs.[[7|, [| which were obtained by the 
SUSY methods and explicit form of zero-energy wave functions. This verifies that 
the SUSY methods give correct results. 

This paper is organized as follows. In Sect. 2, field equations of random-mass 
Dirac fermion is given. In Sect. 3, we shall consider field equations in telegraphic 
configurations of random mass and reformulate the system by introducing the TM 
formalism. In the TM formalism, the eigenvalue problem of the Dirac field reduces 
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to a simple matrix equation. Wave functions are also easily obtained by numerical 
calculation in the TM formalism. Results of numerical calculation are given in Sect. 4. 
We identify localized and extended states from the wave functions. In order to see this, 
our previous studies on non-locally correlated disorder are quite useful. Especially 
it is shown that the DOS obtained by the numerical methods is in good agreement 
with the previous analytical calculation. We also obtain the ensemble averaged Green's 
functions at vanishing energy or for nodeless states and compare them with analytic 
expression in Ref. They are in good agreement though the exact zero-energy wave 
functions used in the analytic calculation are not always normalizable contrary to 
those in the numerical calculation. Section 5 is devoted for discussion. We shall 
discuss relationship between the localization length and the number of states with 
energy below E. Some application of the results for spin chains is also discussed. 



2 Field equations 

We shall consider a Dirac fermion in one spatial dimension with coordinate-dependent 
mass whose Hamiltonian is given by, 

H = J dx^hifj, (2.1) 
h = -ia z d x + m(x)a y , (2.2) 

where a are the Pauli matrices. Defining each component of if) as if) = (u,v), we can 
write the Dirac equation as, 

f — + mix) ) u(x) = Ev(x), 
\ ax J 

( + m(a:) v(x) = Eu{x). (2.3) 

From Eqs. (|2.3|) , we obtain Schrodinger equations, 

f — j— ^ — m'(x) + m 2 (x) j u{x) = E 2 u(x), 



f-^a + m '(x) +m 2 (x)^j v(x) = E 2 v(x), (2.4) 

where prime denotes derivative with respect to the space coordinate. The exact 
solutions to (|2.3|) for E = are easily found though they are not always normalizable, 



u, v oc exp 



± / dx mix) 



(2.5) 



However, solutions for positive eigen-values of E cannot be represented in a simple 
form as in ( |2.5|) . 

On the other hand, when m{x) has a soliton-like configuration, such as[| 

m(x) = mtanh(Amx), (2-6) 

the lowest-energy eigenstate for A = 1 is a "bound" state and explicitly given by 

u(x) = — i— , v(x) = 0. (2.7) 
cosh mx 



We call the configuration (|2.6| ) for A > anti-soliton and that for A < soliton. 
It is one of our conjecture that localized states in the random potential, V r (x) = 
±m'(x) + m 2 (x), are essentially given by linear combinations of the above bound 
states in soliton and anti-soliton configuration of m(x). To see this, we approximate 
the random potential to be a repetition of anti-soliton-soliton like configurations (|2.6|) 



and further, deform it to a series of the step functions by taking the limit A — > oo; 

m (x) = rh(6{x - on) - 1) + ^ fh{6(-x + fa) - 1) + m , (2.8) 

i 3 

where «j and fa are positions of solitons and anti-solitons, respectively (see Fig.l). 
In the subsequent sections, we shall solve the Schrodinger equations in the random 
background fl2.8|). 

1 This configuration of m{x) is a soliton solution in the X(j) 4 Higgs potential of double-well formpO[ . 
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Figure 1: An example of configurations of solitons and anti-solitons, or "rectangular 
barriers" (Eq.P)). 

3 Eigenfunctions 

For nonzero E, solutions of the Dirac equation ( |2.3| ) cannot be expressed in a general 
form as in ( |2.5|) . However, if we take the limit A — » oo, m(x) takes alternatively two 
values m + fh and m — m on intervals whose lengths are positive random variables. 
This drastically simplifies the problem and we can calculate the energy eigenfunctions 
for nonzero E systematically using transfer matrices. In the first half of this section, 
solution in the background m(x) with one step function is considered as a special 
case, and in the second half, it is shown that by connecting so obtained solutions, one 
can obtain solutions for arbitrary patterns of anti-soliton-soliton pairs. 

3.1 Step- functional background 

We study the Dirac equation (|2.3|) with a mass configuration 
wall" at x — 0, 

m(x) = m(29(—x) — 1) + m , 
where 9(x) is the Heaviside function, 

f (x < 0) 
6(x) = { 

[ 1 (0 < x) 
6 



which has a "domain 



(3.2) 



and it is assumed that 

m > mo > 0. 

Substituting ( |3~T|) into (|2.4j) , we obtain a Schrodinger equation with a potential which 



is a combination of the delta function and the step function, 

V(x) = =pm'(x) + m 2 (x) = ±2rhS(x) + ArriQm6(—x) + (m — mo) 2 . (3-3) 

For such a field equation with a delta-function-type potential, the continuity of the 
wave function u(x) leads automatically to the conditions on its derivative u'(x). For 
example, if we rewrite the Schrodinger equation as 

— — - + 2m 5(x) + 4m m 8(—x) u(x) = (E 2 — (m — m ) 2 ) u(x) (3.4) 
y ax 2 J 

and integrate Eq.( |3.4|) with respect to x in the range — e < x < e, we have 

— w'(0 + e) +u'(0 — e) + 2m u(0) + 4m Q rh f dx u{x) = {E 2 — (m — m ) 2 ) / dxu(x). 



X 



(3.5) 

However, the integrations in Eq.( |3.5| ) vanish in the limit e — > due to the continuity 
of u(x). A similar argument applies to v (x) and we have conditions, 

u'(0 + 0) - u'{0 - 0) = 2mw(0), 

v'{0 + 0) - v'{0 - 0) = -2mv(0). (3.6) 

We look for solutions under these conditions. 
Let the energy eigenvalue E satisfy 

m-m > E >0, (3.7) 

then u(x) takes the form 

(3.8) 
(3.9) 

(m-m ) 2 -£ 2 . (3.10) 




From the continuity of u(x), 



A+B=C+D 



and Eg. ( p.6|) gives 
Hence, 

and for x > 0, 

u(x) - 



k 2 {C — D) — K\(A -B) = 2rh{A + B) 



u{x) = Ae Klx + Be~ KlX , for x < 



2m - kx + k 2 2m + Kx + k 2 

-A H B 



(3.11) 



(3.12) 



(3.13) 



2k 2 



*2 



+ 



2m + Kx + k 2 -2m - Kx + k 2 

— A H B 



2k 2 k 2 
If we let m = 0, Ki = k 2 and Eq.( p.l4|) simplifies to 

A e KX + B e~ KX , 



(3.14) 



u(x) = < 



m A m + k 
—A H B 

K K 



e KX + 



-m + k ^ m 

K K 



(x < 0) 
(0 < x) 



(3.15) 



Values of A, B and E are determined by boundary conditions, e.g., u(—L/2) = 
u(L/2) = for a given system size L. 

Our assumption that u(x) remains continuous in the limit A —>■ oo is justified as 
follows. The Schrodinger equation with the potential 



V(x) = ±m'(x) + m 2 (x), 



mix) 



mtanh(Ama;) 



(3.16) 



is transformed to hyper-geometric differential equation by the change of variables, 

1 



1 + e 



-2\mx ' 



(3.17) 



and solved exactly. If we let A tend to infinity, the solution coincides with (|3.14|) . 
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3.2 Transfer-matrix formalism 



In this subsection we restrict m to be zero. Since the wave function u(x) is expressed 
everywhere as 

u{x) = A e KX + B e~ KX , (3.18) 

we represent the eigenfunction in terms of coefficients A and B in what follows. By 
using matrix representation and from ( p. 14 ), the conditions (|3.6|) gives the following 
relation between the coefficients A and B, 



' A(x > 0) 
v B(x > 0) 



T 



A(x < 0) 



y B{x < 0) 



/ m m \ 
1 1 + — — ^ 

K K 



\ 



m m 

K K 



(3.19) 



For soliton instead of anti-soliton, one should replace m with (— m) in Eq. (|3.19|) and 
one obtains 



/ 



A(x > 0) 



-i 



Let us define 



^ B{x > 0) 



R(k, a) = 



T 



( ^na 



A(x < 0) 
B(x < 0) 



(3.20) 





e~ Ka 



K 



(3.21) 



and TM for the configuration of an anti-soliton and a soliton is given as 

R{b)T- l R(a)T 

(1 - 2 ) e K ( a+fc ) + (j) 2 e^ b -^ 0(1 - 



( 
v 



g«(a+b) a) 



gft(a— b) g— fc(a+6) 



kp.\ p -K{a+b) _|_ J ) 2 e «(a-fe) 



= T(a,b) 



( gKa 



2(f) 2 sinh na 2<p( 1 — 0) sinh Ka 



y 20(1 + 0) sinh na e Ka + 20 2 sinh 



(3.22) 



where a is the distance between the soliton and anti-soliton. We impose boundary 
condition on the wave function such that the wave function vanishes as \x\ tends to 
large. This means that the solution, which is proportional to e KX for x — > — oo, does 
not have the factor of e KX for x — > oo. Hence k needs to satisfy 

( 1 ) T(a, b) ( 1 j = e Kb (e Ka - 20 2 sinh «a) = 0. (3.23) 



Since e Kb ^ 0, this means 



K 



2 



= 1 - e - iKa . (3.24) 



?71 



The energy eigenvalue tends to vanish in the limit k — > m, a — > oo, as expected. If 
there are two pairs of anti-soliton and soliton, the boundary condition leads to 

1 ^)T(c,d)T(a,b) ( * j =0, (3.25) 
or 

= e K{a+b+c) + 20 2 [ sinh na sinh nc e~ Kb - e K(b+c) sinh na - e K{a+b) sinh nc ] 

+80 4 sinh na sinh nb sinh kc. (3.26) 

The above argument can be generalized to an arbitrary numbers of anti-soliton-soliton 
pairs readily and we have the following eigenvalue equation, 

/ 



1 )T(e,f)---T(c,d)T(a,b) 



0. (3.27) 



V 



It is easily seen that variables k, m and distances between solitons a, b, c, • • • appear 
always in the combinations rh/n, na,etc. Hence the eigenvalue equation Q3.27 ) is 
invariant under a transformation 

E — > aE, fa — > am, a^a/a (3.28) 

with an arbitrary positive constant a. This property becomes important when we 
discuss relation between the localization length and the energy dispersion as we shall 
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see in later sections. We are not able to solve ( |3.27| ) analytically for arbitrary config- 
uration of pairs of solitons. However it is not so difficult to solve (|3.271 ) numerically. 
We can also easily obtain eigenf unctions after having eigenvalues by using the TM. 
In subsequent sections, we shall give solutions to Eq. (|3.27| ) obtained numerically for 
various multi-soliton configurations m(x). Then we shall compare the results with the 
analytical calculation by the SUSY methods in our previous paper 0. It is straight- 
forward to extend the above formalism for nonvanishing mo- Numerical solutions will 
be given also for m 7^ 0. 

4 Results of numerical calculation 

4.1 Overview 

We solved ( |3.27| ) numerically in various multi-soliton backgrounds, and obtained ex- 
plicit form of corresponding eigenfunctions. We are interested in the relationship 
between randomness and localization of eigenfunctions. Then we compare wave func- 
tions of the eigenstates in quasi-periodical backgrounds and random backgrounds. 

1. Quasi-Periodic Background 

First, we consider a quasi-periodic background in a system of finite size L, 
with almost equal distances between each successive soliton and anti-soliton. 
Eigenfunctions obtained by numerical calculations show twofold oscillations, 
i.e., rapid oscillations occurring at positions of soliton and anti-soliton, and 
slowly oscillating envelope (see the graphs at the top of Fig. 2). The peaks of 
rapid oscillation can be considered essentially the same as the bound state ( pT7|) 
in a single soliton background, because in a background of sufficiently separated 
solitons, eigenfunctions can be approximated by a linear combination of bound 
states located at solitons and anti-solitons. 

The envelope of eigenfunctions is reminiscent of sine curves and the nodes of 
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eigenfunctions coincide with those of corresponding sine curves. Frequencies of 
the envelopes are larger for higher energy eigenvalue and this supports "one- 
dimensional node counting theorem" 0, which states that the number of nodes 
of eigenfunction increases with eigenenergy. 

In quasi-periodic background, all eigenstates extend over the whole system and 
there are no localized eigenstates. This is reminiscent of Bloch's theorem al- 
though we do not impose the periodic boundary condition on the wave functions. 

2. Randomly Distributed Barriers: m = 

We now vary the distances between soliton and anti-soliton according to the 
gaussian distribution, i.e., lengths of rectangular barriers I are subject to the 
gaussian distribution, 



m is set to zero here. In this case, the periodicity of eigenfunction is lost as 
randomness or a in Eq . (f4. 1| ) increases and the envelope has a large peak. This 
means all the states tend to "localize" (see Fig.2). The state of the lowest energy 
E ~ extends over the whole system as it can be expected from the formula 
( |2.5|) . Low-energy states are still extended even at a = 0.4. 

The degree of randomness is controlled by the variance of gaussian distribution 
cr. For larger a, shapes of wavefunctions become sharper and the states become 
more localized. We took a close look at various numerical results and found that 
some of the states have more than one dominant peak in two or three intervals 
between nodes and others have only single "large" peak and are localized in 
an interval between adjacent nodes. Let N(E) denotes the number of states 
from the lowest to the energy E per unit length. Since the width of the interval 
is approximately proportional to l/N(E), we are led to a conjecture that the 




1 



exp{-(7 - fi) 2 /2a 2 }. 



(4.1) 
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n= 1 sigma= 




n= 3 sigma= 




n= 5 sigma= 




n= 1 sigma= . 2 





n= 3 sigma= . 2 





n= 5 sigma= . 2 




n= 1 sigma= . 4 




n= 3 sigma= . 4 




n= 5 sigma= . 4 




n= 1 sigma= . 6 




n= 3 sigma= . 6 




n= 5 sigma= . 6 




Figure 2: The shapes of u(x) in rectangular barriers distributed according to gaussian 
distribution with the average lengths (I) = 2 and fh — 1. In this figure n denotes the 
state of n-th energy from the lowest and m is set to zero. The graphs at the top are 
those in a periodic potential, where the standard deviation a — 0. Here a increases 
from zero at the top to 0.6 at the bottom by 0.2. Eigenfunctions v(x) have similar 
behavours to u(x). 



13 




X 



X 



Figure 3: The shapes of u(x) in rectangular barriers whose lengths I are distributed 
according to exponential distribution with A = 0.25 (Eq. Q4.3|) ). In this figure n 
denotes the state of n-th energy from the lowest, rh and m are set as m = and 
rh = 1. 
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localization length £(E) satisfied 



A 



(4.2) 



N(E) ' 



The proportional coefficient A seems to depend on the distribution of the length 
of rectangular barriers I, m and L, but not on E. 

In Fig. 3, we show eigenfunctions in random rectangular barriers whose lengths 
I are distributed according to the following exponential distribution instead of 



From these results, we can see that the localization around one peak is enhanced 
in higher-energy states. We also verified that as A increases extended states are 
hindered, just as in the gaussian distribution. 

3. For m ^ 

In the above we dealt only with the case of vanishing m®. Here, we take mo ^ 0. 
For quasi-periodic backgrounds, states are still periodic and extended. "Bloch's 
theorem" is valid as in the case of mo = 0. If we vary the length of barrier / 
randomly, states are more localized than in the case of vanishing mo (see Fig. 4). 
Thus we conclude that nonzero m is not necessary for localization, but at least 
it seems to enhance the localization. 

4.2 Comparison with the analytical results: Density of states 

Although Dirac fermion systems with non-locally correlated mass m(x) have sel- 
dom been studied, the case of exponentially distributed m(x) was studied recently 
2 As is well known, there are two different localization lengths in the present system, i.e., the typical 
and disorder- averaged localization lengths (Griffiths phase). In this discussion, £(E) is considerded 
as the averaged localization length, which is defined by the long-distance behaviour of disorder- 
averaged Green's functions, i.e., we expect that the Green's functions are dominated by the states 
with more than one peak||. 




(4.3) 
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X 



X 



Figure 4: The shapes of u(x) in rectangular barriers whose lengths I are distributed 
according to exponential distribution with A = 0.25 (Eq. Q4.3|) ). In this figure n 
denotes the state of n-th energy from the lowest. The parameters m and tuq are set 
as mo = 0.3 and m = 1. 
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by Comtet et.al.H via stochastic method and more systematically by Ichinose and 
Kimura via SUSY methods. Analytical expression of the DOS was obtained there. 
We calculate the DOS numerically by using the TM method and compare the result 
with analytical expression. 

For locally distributed m(x), the density of states with energy E, p{E) = dN X , 
is given by 

as E tends to small. For non-locally correlated random mass m(x)p|], 

1 



Px(E) 



2f llnf |3 



-4g\ 



2f |lnf | 4 



30 f|lnf|3 2#|ln#| 4 #|ln#| 5 

4g I 4p I 4s I 4s I 4s I 4s I , 

4s I 4s I 

where g and A are defined a^ 

[ m(x) m(y) ] cns = £ exp -y|/A). (4.6) 
A 

Parameter A is the correlation length of the disorder. As A — > 0, m(x) and m(y) for 
x y become uncorrelated, the white- noise limit. 

Distribution corresponding to Q4.6|) can be realized by rectangular barriers whose 
height is 

m = Ji, (4.7) 



and width between solitons, I, is distributed asfl 



P{1) = 2l 6XP 



2A 



(U 



3 Definitions of g in Refs.Q and are in fact different by factor 2. In this paper we follow that 

i- 



17 



The result of numerical calculation of 300 pairs of a soliton and anti-soliton is 
shown in Fig.5. We obtained the number of states N num (E) by numerical calculation, 
and we define p num as the averaged density of states around the energy level E as 

PnUE) = jV "--( E + ¥) A£ ^(g-¥) . (4J) 

In order to compare the numerical result of the density of states with energy E with 
the analytical expression p\{E) in (|4~5|), we see the ratio r(E) = p P n ^ E) ■ If the energy 
dependence of these density of states coincides, this ratio r(E) should be constant. 
In Fig.5. r(E) is shown. From this, we can conclude that the energy spectrum of the 
states obtained by the numerical calculation is in good agreement with the analytical 
result ( }4.5| ). Especially, the higher-order expression of p\(E) in ( |4.5| ) gives the better 
agreement. This result also indicates that 300 pairs of a soliton and anti-soliton is 
large enough for the investigation of the system in the low-energy region. 

4.3 Green's functions at vanishing energy 

In this section, we show numerical result of the ensemble averaged Green's functions 
at vanishing energy and compare them with analytical expression in Refs.p, |J. 

The zero-energy wave functions are given by Eq. (|2.5| ). It is not so difficult to 
calculate the ensemble-averaged correlation of them 

W q (x,L) = [\^(x)^(0)\% ns (4.10) 

in the white-noise limit with the system size L||, Hj. There the normalization of 
the wave functions plays an important role because only for specific m(x) they are 
genuine normalizable functions. More explicitly, in solving the Dirac equation we 
imposed the boundary condition such that the wave function decays exponentially 
outside of the random potential. In the analytical calculation in || |J, such kind of 
boundary condition is not imposed on the wave function. 
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Rh o_S U S Y / Rh o_Num 




Order 3 



I — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — Energy 
0.05 0.1 0.15 0.2 0.25 

Figure 5: Energy spectrum of fermion: A comparison between the result of numerical 
calculation and p\{E) of the zero-th order and up to the third order in A, Eq. (|4.5| ). 
Here we set A = 1/6, g = 1, AE = 0.02 and numerical result is averaged over 4000 
trials. The ratio r(E) is shown in this figure. If the numerical calculation agrees with 
Px(E) which is obtained by SUSY methods, the ratio should be constant. The dashed 
line is expected exact result. 
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In general it is expected that the above correlation functions exhibit complex 
multi-fractal scaling as in the qunatum Hall state, 

W q (x,L) ~ L- d - T ^\x\- y{q) . (4.11) 
For the present system, W q (x, L) is obtained as follows for L — > oo, 

W,(,,L)~*<£> (4.12) 



ty(l)= Ptft , r,,„ e- rf . (4.13) 



with 

y (TTp) 4 ' 

For large x 

W q (x,L) ~ l/(x 3/2 L) (4.14) 

and therefore it is expected that r(g) = and y(q) = 3/2. 

We numerically calculated the correlation functions W q (x,L) with various values 
of q. To this end, we focused on the states of nodeless wave function. The results 
are shown in Figs. 6 and 7 for q = 0.5, 1.0, 1.2 and 1.5 as well as the analytical 
result Eqs.( f4.12|) and ( 14.131) . It is obvious that the numerical calculations are in 



good agreement with the analytical calculations. Especially we can conclude that the 
critical exponent of the multi- fractal scaling is obtained as y(q) =3/2 independently 
of q. (It seems that y(q) is smaller than 3/2 for small q. But we expect that y(q) 
with any q approaches to 3/2 for larger systems.) 

We also study L dependence of W q (x,L). The result is shown in Fig. 8. Accord- 
ing to this, it seems that W q (x,L) is not inversely proportional to L and moreover 
r(q) 7^ 0, contrary to the analytical calculation by Balents and Fisher 0. We further 
calculated W q (x,L) for larger system size L and smaller g, and obtained similar re- 
sults. However it is still possible that the system size L is not large enough in our 
numerical calculations. Therefore the above is not a definite conclusion and more 
systematic calculation is desired. (In the analytical calculation 0, L is assumed to be 
large enough.) 
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We should again remark here that Eqs.( |4.l2| ) and ( |4.13| ) are obtained from ( |2.5| ) 



some of which are not normalizable. On the contrary, the wave functions in the 
numerical calculation are all normalizable. Then strictly speaking, one cannot expect 
coincidence of the analytical and numerical results. 

5 Discussion 

From the results of the DOS and the Green's functions, we conclude that 300 pairs 
of a soliton and anti-soliton can be regarded as a good approximation for infinite-size 
system. Hence it is legitimate to discuss behaviors of the localized states based on 
properties of numerically obtained solutions. 

Results obtained from the numerical calculations can be summarized as follows. 

• For a quasi-periodic background, "Bloch's theorem" is valid and all the states 
extend over the whole system, regardless of the value of m . The envelopes of 
wave functions are similar to sine curves. 

• For m = 0, on varying the length of rectangular barriers randomly, states 
begin to localize, within two or three intervals between nodes. Some of them, 
especially the lowest-energy states, has more than one peak which separate with 
each other rather long distance. 

• For nonzero value of mo, the states tend to localize in narrower intervals than 
in the case of mo = 0. Hence nonzero mo enhances, but is not necessary for the 
localization. 

We conjectured that the averaged localization length £,(E n ) for the n-th state from 
the lowest satisfies 



since states in randomly distributed rectangular barriers are localized in two or three 
intervals between nodes and length of each interval 1(En) is proportional to 1/N. 
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Figure 6: Multifractal behaviour of bound-state wavefunction (1) : The x dependence 
of W q (x, L) based on numerical and analytical calculations is shown. Difference be- 
tween two results is also shown. To calculate W q (x, L), we used the u-component 
of the doublet ip = (u,v), namely W q (x,L) num = [\u(x + l)u(l)\] ens . To obtain 
W q (x, L) num numerically, we averaged W q over 2800 configurations of the random- 
varying mass. We also took an average with respect to / in [\u(x + l)u(l)\] ens . 
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Figure 7: Multifractal behaviour of bound-state wavefunction (2) : The method 
to obtain W q is the same as we mentioned in the caption of Fig.6. We can see 
that the inclination of W q with q = 1.2 and 1.5 is almost 3/2 in the region around 
log(x/L) ~ —1. (We set L(the system size)= 40, A = 1/512 and g = 1/2 (almost 
white- noise case) in the numerical calculation, and we set L = 40 and g — 1/2 in the 
analytical calculation.) 
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Figure 8: Multifractal behaviour of bound-state wavefunction (3) : L (the system size) 
dependence of W q (x,L) for q = 0.5 and 1.0 is shown. Here we set A = 1/512 and 
g — 1/2, and W q is averaged over 1600 backgrounds of the random mass. To calculate 
W q , we used only the u component, though we verfied that the v component behaves 
similarly. 
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Here, A in depends on the distribution of lengths of rectangular barriers I, fh 
and L, but not on energy E. One of supports for this conjecture is given by the exact 
result for locally correlated random mass m(x), since the number of states below 
energy E is given as 

N(E) oc (5.2) 

and localization length is given as 

£(£)cx(ln£) 2 , (5.3) 

hence ( |5.1| ) is satisfied in this case. (Note that N(E) and n are related as N = 
n/L.) Very recently, we obtained the localization length of the ensemble- averaged 
one-particle Green's functions by using the SUSY methods^, 

£(£) = ^ + Yg\) + ((^) 2 )- ( 5 - 4 ) 

From this result and the DOS in Eq. (f4.5| ), we can verify 



£{E)N\{E) oc 1 + 0((g\) 2 ). (5.5) 

We also performed numerical calculations with the TM method under Dirichlet 
boundary condition, which is u(0) = u(L) = 0. However qualitatively features of the 
localization are the same. Another important insight into the Dirac fermion system 
with random mass m(x) obtained by the present study is that the n-th state from 
the lowest is related to the lowest-energy state in system of size 1/N(= L/n), since a 
state localized inside an interval of length 1/N must be realized as the ground state 
of the system of length 1/N. 

The ensemble-averaged Green's functions were also calculated. They are in good 
agreement with the analytical results. Especially, though each wave function exhibits 
rapid oscillations, the ensemble-averaged Green's functions have smooth behaviour. 

Finally let us comment on the spin systems closely related to the random-mass 
Dirac fermions. As discussed in Refs.||, §j low-energy excitations in doped spin-ladder 
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or spin-Peierls systems are described by the random-mass Dirac fermions. Undoped 
cases of the above systems have an energy gap, and by doping there appear midgap 
states [p]]. Results obtained in this paper suggest that properties of the midgap states 
strongly depend on the randomness of the doping, i.e., almost all midgap states 
are extended for quasi-periodic doping, whereas states tend to localize for random 
doping. In terms of spin-Peierls system (like the compound Cui-^Zn^GeOs [0), 
the localization length £(E) is related with the correlation length of spins, since the 
existence or absence of our fermion is related to the up or down of the spins via 
Jordan- Wigner transformation Using our results one will be able to find that 
randomly doped spin-Peierls compounds have the short correlation length of spin. It 
is very interesting to see the above behavior by experiments in which the distribution 
of impurities is controled. We are now calculating Green's functions in having more 
definte prediction for the above spin systems. The result will be reported in a future 
publication ||. 
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